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This article reports an automated approach to the evaluation of higher-order terms of QED perturbation to 
anomalous magnetic moments of charged leptons by numerical means. We apply this approach to tenth-order 
correction due to a particular subcoUection of Feynman diagrams, which have no virtual lepton loops. This 
set of diagrams is distinctive in that it grows factorially in number as the order increases, and also each of the 
diagrams holds quite a large number of subtraction terms to be treated along renormalization procedure. Thus 
some automated scheme has long been required to evaluate correctly this class of diagrams. We developed a fast 
algorithm and an implementation which automates necessary steps to generate from the representation of each 
Feynman diagram the FORTRAN codes for numerical integration. Currently those diagrams of tenth order are 
being evaluated. 



1. Introduction 

The anomalous magnetic moment of electron, 
also called the electron g — 2, is one of the fun- 
damental quantities of particle physics. Since its 
discovery in 1947 (jj it has been measured with 
steadily increasing precision |2I3| . The best val- 
ues of g—2 of the electron and the positron avail- 
able in the literature ^ 

ae- (exp) = 1 159 652 188.4 (4.3) x lO"^^ , 
ae+(exp) = 1 159 652 187.9 (4.3) x 10"^^ 

were obtained by the Penning trap experiment. 

The electron g—2 (denoted by = ^{g — 2)) is 
explained almost entirely by the electromagnetic 
interaction between electron and photon alone, 
and thus it has provided the most stringent test 
of QED. An important by-product of the study of 
the electron g~2 is that currently the most precise 
estimation of the fine structure constant a can be 
obtained by combining the measurement and the 
theory of Og, which yields 0] 

a-i(ae) = 137.035 998 834 (12)(31)(502), 
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where the uncertainties 12 and 31 are due to 
and terms, and 502 comes from the experiment 

o . 

At present a new experiment on the electron 
g — 2 is being carried out by a Harvard group, 
which will reduce the measurement uncertainty 
substantially 5 . It will enable us to test the va- 
lidity of QED to a very high degree and to deter- 
mine a to an unprecedented precision of 7 x 10~^° 
or better. Of course such a feat requires availabil- 
ity of the theoretical calculation of matching pre- 
cision. The present issue is that the uncertainty 
due to term may become major source of sys- 
tematic errors. Therefore a reliable theoretical 
estimation of corrections is urgently required. 

From the viewpoint of obtaining the correc- 
tions the numerical integration approach is the 
only practical choice at present. We employ the 
numerical evaluation scheme developed by one of 
the authors (T. K.) and Cvitanovic |6I7| . 

The contribution to the term of the elec- 
tron g—2 comes from 12672 vertex- type Feynman 
diagrams, which can be categorized into 6 sets 
according to their structures and classified fur- 
ther into 32 gauge- invariant sets .?< !) . None of 32 
sets is dominant so that all must be evaluated. A 
particularly difficult one is Set V, a huge set con- 
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Figure 1. A typical diagram of Set V, which is 
represented as {(0, 2)(1, 4)(3, 6)(5, 8)(7, 9)}. 



sisting of 6354 vertex diagrams, aU of which have 
radiative corrections only due to virtual photons. 
We begin by exploiting the equation derived from 
the Ward-Takahashi identity, 



9A^(p, q) 



which relates the sum of 9 vertex diagrams to one 
self-energy-like diagram. This relation together 
with the time-reversal symmetry of QED enables 
us to reduce the number of independent integrals 
of Set V diagrams drastically from 6354 to 389. 

The difficulty of Set V also stems from the 
fact that many of them have very large num- 
ber of UV and IR divergences, e.g. maximally 
47 UV counterterms are required for some dia- 
grams and many of the others require more than 
20 UV counterterms. This makes the previous 
approach highly impractical since it runs into an 
extremely sever logistic problem. Toward this dif- 
ficult problem we present our solution by an au- 
tomated scheme for code generation, which en- 
ables us to obtain renormalized integrals for all 
diagrams of Set V that are ready to be integrated 
by numerical means 

2. Diagrams without lepton loops 

Set V diagrams summed by W-T identity are 
one-particle-irreducible (IPI) self-energy-like dia- 
grams of 2nth order which have no closed lepton 
loops. A typical diagram is shown in Figure ^ 
They are given by a path consisting of all lep- 
ton lines emanating from the incoming lepton and 
terminating at the outgoing lepton. 2n interac- 
tion vertices lie on the path, which are labeled by 
Wo , wi , . . . , V2n-i from the left to the right in our 
conventions. A photon line labeled by k connects 
the vertices Vi^ and Vj^ at its ends. Therefore 



a diagram is uniquely specified by a single-line 
expression {(vi^,Vj^)^...} which represents the 
pattern of pairings how photon lines connect the 
vertices if we adopt a proper ordering of pairs. 

All Set V diagrams are generated systemati- 
cally by finding every possible pattern of pairings 
that satisfies IPI condition, which reduces to a 
matter of simple combinatorics. Due to their sim- 
ple structure the graph-theoretical notions, espe- 
cially the independent set of closed paths on the 
diagram, are easily identified. It is crucial in con- 
structing the integrand. Another significant fea- 
ture is that UV-divergent subdiagrams are eas- 
ily found as single segments on the lepton path. 
The UV-divergent structure of the diagram will 
be identified exactly by forests constructed by the 
combinations of the subdiagrams. Those proper- 
ties are quite efficiently handled in an automated 
manner. 

3. Numerical integration formalism 

The anomalous magnetic moment is given 
by the static limit of the magnetic form factor 
^2(9^)1^2^0- ^^'^ contribution from a 2nth order 
Feynman diagram Q is expressed by an integral 
over n loop momenta {kr}: 



-Y.g = {ie) 
I 



n 



n 



where Or^i — r^5^ is a diagram-specific product 
specified by the pairing pattern in the case of Set 
V diagrams. The momentum integration is ex- 
actly carried out and the amplitude is converted 
into an integral over 13-dimcnsional Feynman pa- 
rameter space: 



Vtt/ 4" ^ 



n-1 



Fo [Bij , , 



Fi{Bij,Aj,Cij] 

U'iyn-2 



The result is expressed symbolically as a func- 
tion of quantities [/, Bij , Aj , and V (called as 
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building blocks), which are homogeneous polyno- 
mials of Feynman parameters. The explicit forms 
of U and Bij are determined by the underlying 
topological structure of the Feynman diagram Q 
[TT] . Aj , Cij and V also reflect the structure of 
the diagram and are obtained from U and Bij. 

The amplitude constructed above is divergent 
in general, and the divergences must be removed 
before carrying out the integration numerically. 
We adopt the subtractive on-shell renormaliza- 
tion. In order to perform renormalization numer- 
ically we employ a strategy in which we prepare 
the subtraction terms as integrals over the same 
domain of the original unrenormalized amplitude, 
and carry out point-wise subtraction so that the 
singularities of the original integrand are canceled 
point-by-point on the Feynman parameter space. 

To achieve this we adopt the following interme- 
diate renormalization scheme. We first subtract 
the UV-divergent part of the unrenormalized in- 
tegrand, which is identified by simple power- 
counting rules. This procedure is formulated as 
JT-operation. By construction the subtraction 
term factorizes exactly into UV-divergent part of 
the renormalization constant Lg associated to the 
UV-divergent subdiagram S and the lower-order 
g — 2 term Mg/g. The factorization property is 
significant in the present formulation in that it 
allows successive application of JT-operations. 

The whole subtraction terms are given by Zim- 
mermann's forest formula jj^l- Each source of 
UV-divergence is related to a forest, a set of non- 
overlapping UV-divergent subdiagrams. The sub- 
traction term associated to a forest is constructed 
by applying JT-operations to the unrenormalized 
integrand successively for the subdiagrams in a 
proper order. The advantage of the forest ap- 
proach is that it is readily translatable into code 
generation. The forests are given by the combi- 
nation of the subdiagrams, so the complete iden- 
tification of UV-divergent parts are obtained by 
purely combinatorial procedure. Thus it is quite 
efficiently implemented in terms of forests, which 
enables us to obtain fully UV-renormalized am- 
plitude of the Feynman diagram Q. 

The diagram may have IR-divergences when it 
contains self-energy subdiagrams. The subtrac- 
tion of IR divergences is achieved by /-operation 



in a similar manner as that of UV divergences. 

Counterterms thus constructed can be identi- 
fied with only the divergent part of the renormal- 
ization constants so that the result of the above 
step is not fully equivalent to the standard on- 
shell renormalization. To complete the calcula- 
tion the difference between the full renormaliza- 
tion and the intermediate renormalization must 
therefore be evaluated by summing up all sub- 
traction terms. This step is called the residual 
renormalization, which is also much complicated 
for the higher-order calculations and is to be au- 
tomated in the future work. 

4. Automated flow of calculation 

The flow of the process to generate the numer- 
ical integration code for evaluating an individual 
diagram is schematically depicted as Figure El 

The information of a diagram is provided in a 
single-line representation indicated by the rectan- 
gular box at the upper-left corner. It enables us 
to construct the unrenormalized amplitude deter- 
mined directly by the pairing pattern of the dia- 
gram, which is converted into the form of Feyn- 
man parametric integrals in terms of building 
blocks by analytically performing the integration 
over the loop momenta with the help of algebraic 
manipulation program, FORM 

The explicit form of the build blocks are de- 
termined from the underlying topological struc- 
ture. They are obtained by referring to the funda- 
mental set of closed paths on the diagram, which 
can be extracted from the diagram representa- 
tion. This step also incorporates algebraic ma- 
nipulations performed by MAPLE. 

The UV divergences of the diagram are iden- 
tified as Zimmermann's forests, which are ob- 
tained by the non-overlapping combinations of 
UV-divergent subdiagrams once the complete set 
of those subdiagrams are found from the di- 
agram representation. The subtraction opera- 
tion is achieved by if-operation, in which simple 
power-counting rules are applied to the original 
unrenormalized integrand and the building blocks 
to result in the proper form of subtraction term. 
A subtraction term corresponding to a particular 
forest is constructed by the successive application 
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Figure 2. Flow of process to generate the numerical integration codes from the diagram representation. 



of i(r-operations, each of which is associated to the 
subdiagram in the forest. 

The IR divergences remaining in the individ- 
ual diagram should also be subtracted away by 
/-operations in a similar manner as those of UV 
divergences, though this subject is not covered in 
the present article. 

Finally, the (intermediate-) renormalized am- 
plitude constructed from the original amplitude 
and the set of subtraction terms is turned into 
a FORTRAN code. It is readily processed by 
the numerical integration system such as VEGAS 
an adaptive Monte-Carlo integration rou- 
tine. 

5. Current status and concluding remarks 

In this article we presented an automated 
scheme of code generation for evaluating higher- 



order QED corrections of the electron anomalous 
magnetic moment by numerical means. We con- 
structed an algorithm and concrete procedure to 
obtain UV-renormalized amplitudes for a partic- 
ular set of diagrams without lepton loops in an 
automated manner. It should be noted that the 
scheme itself is applicable to an arbitrary order 
of diagrams of this particular type, though our 
practical concern is to evaluate the tenth-order 
corrections. 

We obtained a program to list up all the topo- 
logically distinct diagrams, and to identify the 
UV divergent structure of each diagram in terms 
of forest structures. We implemented our au- 
tomated procedure as a set of Perl programs 
with the help of symbolic manipulation systems, 
FORM and MAPLE. From a single-hne repre- 
sentation of a diagram it generates numerical in- 
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tegration codes in FORTRAN, which are ready 
to be processed by VEGAS, an adaptive Monte- 
Carlo integration routine. 

The programs have been tested for lower-order 
diagrams and confirmed that they reproduce the 
codes for the sixth-order and eighth-order dia- 
grams previously constructed. They are now be- 
ing applied to tenth-order diagrams. At present, 
the diagrams which have only vertex renormaliza- 
tion were processed and test runs were performed. 
Those diagrams corresponds to 2232 vertex dia- 
grams among 6354 Set V diagrams of tenth-order. 
Crude evaluation showed no sign of divergent be- 
havior, which confirms that our scheme is work- 
ing well as expected. They are currently put to 
production runs. 

Typical size of the numerical integration code 
of a diagram amounts to about 80,000 lines of 
FORTRAN program. It takes 10-20 minutes to 
generate the code on an ordinary PC. The run- 
ning time to evaluate 10^ samplings times 20 iter- 
ations takes about 5-7 hours on 32 CPU cluster. 
It seems that the estimated computational cost 
required to accomplish the whole numerical in- 
tegration in a few percent of uncertainty remain 
within the manageable scale with up-to-date com- 
putational facilities. 

The remaining 4122 diagrams have not only 
UV-divergent self-energy subdiagrams but also 
IR divergences. The simplest way to deal with the 
IR problem is to give a small mass A to photons, 
which requires no further work on the automat- 
ing code, and is being pursued as a first step. To 
obtain a result independent of A it is necessary 
to incorporate IR subtraction terms in a manner 
similar to that of UV counterterms, which is left 
as our future issue. 
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